The vector approach

I will use a vector approach to this assignment, where I will use the various layers to create a set of new features (multipolygons) that represent suitable areas by each criterion. The intersection of all these features is the set of suitable areas for wind development.

Read in data

We’re given three data files for this assignment: basemaps.gpkg, inputs.gpkg, and Wind2002_which_one. Let’s see what’s in inputs.gpkg:

# Create a list of layer names of inputs.gpkg and display them
input_layers <- st_layers("HW3/data_hw3/inputs.gpkg")$name
input_layers
## [1] "Airports"        "CenterLines2008" "FireLRA"         "Parcels"        
## [5] "FireSRA"         "County"          "Urban2010"

We could read these into a single sf object using a for loop, but in this case I’ll keep them as separate objects, so I’ll just read them in one by one. The read_sf function will automatically recognize the format of the file to be read in, so we just need to specify the layer (using the list of layer names created in the previous code chunk).

# Read in input layers individually
airports <- read_sf("HW3/data_hw3/inputs.gpkg", layer = input_layers[1])
roads <- read_sf("HW3/data_hw3/inputs.gpkg", layer = input_layers[2])
fireLRA <- read_sf("HW3/data_hw3/inputs.gpkg", layer = input_layers[3])
parcels <- read_sf("HW3/data_hw3/inputs.gpkg", layer = input_layers[4])
fireSRA <- read_sf("HW3/data_hw3/inputs.gpkg", layer = input_layers[5])
county <- read_sf("HW3/data_hw3/inputs.gpkg", layer = input_layers[6])
urban <- read_sf("HW3/data_hw3/inputs.gpkg", layer = input_layers[7])

The roads layer (along with all the other layers) has the following CRS:

# Check CRS
st_crs(roads)$input
## [1] "NAD83 / California Albers"

Finally, we read in the wind raster data. For this we need to use the raster function from the raster package.

# Read in raster
wind_raster <- raster("HW3/data_hw3/Wind2002.tif")

Wind

One thing to note about the wind data is that it is a categorical variable, with seven wind classes (1, 2, …, 7). We see this from looking at the unique values in the data:

# Unique values of the wind speed data
unique(getValues(wind_raster))
## [1] NA  1  2  3  4  5  6  7

Wind speeds of at least 6.4 m/s correspond to class 3 or above.

I now convert the raster data to polygons. A couple of things to note here. I use the function rasterToPolygons with the option to dissolve, which gives me just one multipolygon per wind speed class. I then use st_as_sf to convert it to an sf object. I’m also going to rename the wind variable.

# Convert raster data to polygons, then turn that into sf object
wind_sf <- wind_raster %>%
  rasterToPolygons(dissolve = TRUE) %>%
  st_as_sf()

# Rename wind speed column
colnames(wind_sf)[1] <- "wind_speed"

It doesn’t hurt to check that this object has the same CRS as the ones previously loaded:

st_crs(wind_sf) == st_crs(roads)
## [1] TRUE

Now we can map this:

# Map
tm_shape(wind_sf) +
  tm_fill(col = "wind_speed", title = "Wind speed") + # color by wind speed class
  tm_shape(county) + # include county boundaries
  tm_borders() +
  tm_legend(position = c("right", "top"), legend.bg.color = "white", legend.bg.alpha = 0.7) +
  tm_layout(main.title = "Wind speeds: vector data", frame = FALSE, main.title.size = 1.5)

Since we require wind speeds above 6.4 m/s, I filter for wind classes 3 through 7.

wind_suitable <- wind_sf[3:7,]

We can map this now.

tm_shape(county) +
  tm_borders(col = "black", lwd = 1.5) +
  tm_shape(wind_suitable) +
  tm_fill(col = "lightgreen") +
  tm_layout(main.title = "Wind speed of at least 6.4 m/s", frame = FALSE, main.title.size = 1.5)

Roads

This requirement states that sites need to be within 7.5 km of a major road. The roads data comes in the form of of MULTILINESTRING (you can check this with st_geometry_type(roads, by_geometry = FALSE)), so we can just make a buffer around those lines. I’ll also clip the roads buffer to the county so it doesn’t unnecessarily extend outside it.

# 7.5 km buffer around roads
road_suitable <- st_union(st_buffer(roads, dist = 7500))

# Clip roads buffer layer to the county
road_suitable <- st_intersection(road_suitable, county)

We can plot this to make sure it looks reasonable.

tm_shape(county) +
  tm_borders(col = "black", lwd = 1.5) +
  tm_shape(road_suitable) +
  tm_fill(col = "lightgreen") +
  tm_shape(roads) +
  tm_lines(col = "red", lwd = 0.5) +
  tm_layout(frame = FALSE) +
  tm_add_legend(type = "line", col = "red", labels = "Road") +
  tm_add_legend(type = "fill", col = "lightgreen", labels = "Within 7.5 km of road")

Airports

We need to identify areas that are more than 7.5 km from an airport. This is equivalent to being inside the county but outside an airport buffer. Now we’ll make a buffer and dissolve it into a single multipolygon.

airport_buffer <- airports %>%
  st_buffer(dist = 7500) %>%
  st_union()

The operation we need to find suitable areas here is st_difference, which, with arguments \(x\) and \(y\) gives us the area that is in \(y\) but not in \(x\).

airports_suitable <- st_difference(county, airport_buffer)
## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries

We can plot this to make sure we got it right:

tm_shape(county) +
  tm_borders(col = "black", lwd = 1.5) +
  tm_shape(airports) +
  tm_polygons(col = "red") +
  tm_shape(airports_suitable) +
  tm_fill(col = "lightgreen") +
  tm_layout(frame = FALSE) +
  tm_add_legend(type = "fill", col = "red", labels = "Airports") +
  tm_add_legend(type = "fill", col = "lightgreen", labels = "> 7.5 km from airport")

Urban area

Similar to the above, we need to identify areas not within 1 mile of an existing urban area. As for airports, we make a buffer and then take the difference between the buffer and the county.

urban_buffer <- st_union(st_buffer(urban, dist = 1609)) # 1609 meters in a mile
urban_suitable <- st_difference(county, urban_buffer)
## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries

And plot it as a sanity check:

tm_shape(county) +
  tm_borders(col = "black", lwd = 1.5) +
  tm_shape(urban_suitable) +
  tm_fill(col = "lightgreen") +
  tm_shape(urban) +
  tm_fill(col = "red") +
  tm_layout(frame = FALSE) +
  tm_add_legend(type = "fill", col = "red", labels = "Urban area") +
  tm_add_legend(type = "fill", col = "lightgreen", labels = "> 1 mile from urban area")

Fire

For Fire LRA, the data are already filtered to only include areas with very high fire risk. For fire SRA, I’ll filter for features with HAZ_CODE = 3. I will take the union of the two, make it a single multipolygon, and take the difference with the county.

# Filter for haz code in Fire SRA
fireSRA_haz3 <- filter(fireSRA, HAZ_CODE == 3)

# Make fire LRA a single multipolygon
fireLRA <- st_union(fireLRA)

# Take difference with county
fireLRA_suitable <- st_difference(county, fireLRA)
## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries
# Make fire SRA a single multipolygon
fireSRA_haz3 <- st_union(fireSRA_haz3)

# Take difference with county
fireSRA_suitable <- st_difference(county, fireSRA_haz3)
## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries

We can now look at the two types of fire zone.

tm_shape(county) +
  tm_borders(col = "black", lwd = 1.5) +
  tm_shape(fireLRA) +
  tm_fill(col = "red") +
  tm_shape(fireSRA_haz3) +
  tm_fill(col = "orange") +
  tm_layout(frame = FALSE) +
  tm_add_legend(type = "fill", col = "red", labels = "Fire LRA zones") +
  tm_add_legend(type = "fill", col = "orange", labels = "Fire SRA zones")

Public land

We want to exclude any areas that are on public land. Public land in the parcels data have a USECODE that is 6000 or higher, so we can first exclude that to create a public lands layer.

public_land <- filter(parcels, USECODE >= 6000)
public_land <- st_union(public_land)

As before, we take the difference in order to find the areas that are not on public land, and make it a single multipolygon.

public_land_suitable <- st_difference(county, public_land)
## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries

Aaand let’s map this.

tm_shape(county) +
  tm_borders(col = "black", lwd = 1.5) +
  tm_shape(public_land) +
  tm_fill(col = "red") +
  tm_shape(public_land_suitable) +
  tm_fill(col = "lightgreen") +
  tm_layout(frame = FALSE) +
  tm_add_legend(type = "fill", col = "red", labels = "Public land") +
  tm_add_legend(type = "fill", col = "lightgreen", labels = "Outside public land")

Final list

Now we have six layers where wind development is suitable by corresponding criteria. To find out where wind is suitable by all criteria, we find the intersection of all these layers.

suitable_all <- st_intersection(wind_suitable, road_suitable, airports_suitable,
                                urban_suitable, fireLRA_suitable, fireSRA_suitable,
                                public_land_suitable)
## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries

Now we have one more thing to do, based on this part of the assignment:

I take this to mean that each suitable polygon has to be at least 4 hectares. What I have in my suitable_all layer is a multipolygon with only one layer, so I’ll transform (“cast”) this into separate polygons.

# Turn suitable areas into individual polygons
suitable_all <- st_cast(suitable_all, "POLYGON")
## Warning in st_cast.sf(suitable_all, "POLYGON"): repeating attributes for all
## sub-geometries for which they may not be constant

Now let’s calculate the area for each of these polygons, and then create a list of polygons that are at least 4 hectares. Since we also want to rank the polygons, I’ll select the polygons with the highest wind speed class (7), then rank by area, and pick the top ten.

# Calculate area of each polygon
suitable_all$area <- drop_units(st_area(suitable_all))

# Top ten
final_list <- suitable_all %>%
  filter(area >= 40000) %>% # at least 4 ha
  filter(wind_speed == 7) %>% # only highest wind speed class
  arrange(desc(area)) %>% # sort by area
  head(10) # select the first 10

Now we can have a closer look at these top 10 areas (in red), with the full extent of suitable areas in green. In this case, I’ll use an interactive map where you can zoom in and out, choose your preferred base layer, etc.

tmap_mode("view")
tm_shape(suitable_all) +
  tm_fill(col = "green", alpha = 0.5) +
  tm_shape(final_list) +
  tm_polygons(col = "red")